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Hydraulic parameters of rock mass are the most effective 
factors that affect rock mass behavioral and mechanical 
analysis. Aforementioned parameters include intensity and 
density of fracture intersections, percolation frequency, 
conductance parameter and mean outflow flowrate which 
flowing perpendicular to the hydraulic gradient direction. In 
order to obtain hydraulic parameters, three-dimensional 
discrete fracture network generator, 3DFAM, was developed. 
However, unfortunately, hydraulic parameters obtaining 
process using conventional discrete fracture network 
calculation is either time consuming and tedious. For this 
reason, in this paper using Artificial Neural Network, a tool 
is designed which precisely and _ accurately estimate 
hydraulic parameters of discrete fracture network. 
Performance of designed optimum artificial neural network 
is evaluated from mean Squared error, errors histogram, and 
the correlation between artificial neural network predicted 
value and with discrete fracture network conventionally 
calculated value. Results indicate that there is the acceptable 
value of mean squared error and also a major part of 
estimated values deviation from the actual value placed in 
acceptable error interval of (-1.17, 0.85). On the other hand, 
excellent correlation of 0.98 exists between the predicted and 
actual value that proves the reliability of the designed 
artificial neural network. 
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1. Introduction 


In designing infrastructures such as dams, tunnels, and caverns which involved in the rock mass, 
hydraulic behavior and connectivity of fractures play important role. The significance of 
aforementioned parameters will be highlighted if the project site placed in zones with high pore 
pressure like underground caverns and tunnels of dams, underground hydrocarbon reservoirs and 
as well as host zones of nuclear waste. The predominantly Hydraulic behavior of discontinuous 
media such as rock mass is controlled by geometry and connectivity of fractures [1-5]. In rock 
mass hydraulic behavior study, due to inaccessibility and direct immeasurability of connectivity 
and flowrate, rock mass fracture network simulating using statistical fracture data such as 
orientation, size and spatial location and density of fracture, which named Discrete fracture 
network(DFN) model, is accepted and common approach [6-9]. Fracture network modeling is an 
important part of the design and development of natural energy resource systems including 
geothermal and petroleum reservoirs and aquifers [10-14]. Robinson studied fracture 
connectivity by percolation theory and numerical methods in 2D and tried to generalize the 
results to the 3D state. He considered fractures as segments which are positioned randomly in 
space, and their orientation follows a distinct distribution function [15]. In discrete fracture 
network model, percolation theory was used by several researchers such as Englman et al., 
Charlaix et al., Long and Witherspon, Stauffer. These researchers used the theory to study effects 
of trace length and size of fractures on the amount of fracture connectivity in 2D [3,16—18]. Bour 
et al. and Bour & Davy, assuming an exponential distribution function, studied different 
properties of fracture connectivity by percolation theory in 2D and 3D [19,20]. Darcel et al., 
investigated fracture connectivity in a discrete fracture network, using fractal correlation 
coefficients, by theoretical and numerical methods [21]. Due to the extreme importance of 
underground cavities for nuclear waste disposal, and to overcome uncertainties of seepage 
situation in the studied area, Xu et al. introduced fracture connectivity index as a parameter to 
quantify the state of fracture connectivity of a region [22]. Xu et al. investigated fracture 
connectivity in geothermal systems by discrete fracture network and studied concept and effect 
of this parameter on the behavior of these systems [23]. During the process of developing a 
discrete fracture network, they used an integrated method of Markov Chain- Monte Carlo. 
During their study on fracture connectivity parameter, Fadakar et al., suggested a new parameter, 
referred to as field fracture connectivity, which represents characteristics of true degree of spatial 
fracture connectivity in geothermal systems [24]. They provided some relations between field 
fracture connectivity, P_21 parameter and density of discontinuities, as well. Discrete fracture 
network is the main element of the study in all of the above methods. However, it should be 
noted that, when using discrete fracture network by generating compact models to achieve a 
representative element volume model, the study of discrete fracture network regarding hydraulics 
becomes very cumbersome. Considering these conditions, obtaining a faster but considerably 
accurate solution is a requirement. 


Addressing this need, it is attempted in this article to design an optimal neural network. Using 
appropriate structural properties such as proper cell model, the exact determination of suitable 
weighted matrix, training algorithm and an activation function, hydraulic parameters intended for 
discrete fracture network were generated receiving input parameters. Application of artificial 
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neural network (ANN) has a long history of use in different aspects of rock mechanics. Most of 
the investigations use the artificial neural network to predict required physical and mechanical 
properties. Ji et al. used ANN as an assistant tool in back-analysis method to calculate 
permeability tensor [25]. Benardos et al. applied ANN to predict TBM! performance [26]. 
Bowden et al. and Kingston et al. used ANN in water resources studies [27,28]. Kung et al. and 
Kim et al. used capabilities of ANN to estimate ground settlement due to underground 
excavations [29,30]. Padmini et al. used an integrated technique of Fuzzy- neural network for 
calculation and estimation of the ultimate bearing capacity of shallow foundations [31]. 
Kulatilake et al. and Menjazi et al. used ANN to predict mean size of rock fragments resulting 
from explosion operations in open mines [32,33]. Nejati et al. and Gokceoglu et al. used ANN to 
obtain deformation modulus of rock mass [34,35]. Aforementioned research works are only part 
of the worked paper which show the applicability of ANN in a different field of rock mechanics. 
This work with several other works indicates the great power of ANN in engineering disciplines. 
For this reasons, ANN is selected as a powerful predicting tool to overcome problems of robust 
calculations of DFN. 


2. Artificial neural network (ANN) 


Inspired by the nervous system of living organisms, Artificial Neural Networks have various 
capabilities. Neural networks are in fact parallel processors that perform a large amount of 
computation on input data. The structure of the neural network is composed of units known as 
neurons and neural layers that offer great flexibility and a large degree of freedom. Based on the 
generation of smart algorithms operation, these networks have a great capacity for correction and 
generalization. Many different models of neural networks were presented with various structures 
and algorithms. However, despite different structures, they are all similar regarding basic 
characteristics such as learning ability and versatility, generalizability, parallel computing, 
robustness, and general approximation capability [36]. 


Every neural network has three structural features in all of which the structure self-organization 
is evident: neural cell model (function type), neural network structure (topology type), and 
learning in the neural network (learning type). The neural cell model is determined based on the 
transfer function. The real output of each neuron depends on the specific transfer function that is 
selected, and it must satisfy the criteria required by the problem that the neural cell is used to 
solve [36]. 


All the elements of the input vector are multiplied by their weight in a neural network and are 

given to the transfer function or the activation function after summation by the bias of each 

neuron. Each layer of the neural network has an input vector, a bias vector, a weight matrix, and 

finally, an output vector which acts as the input vector for the next layer. In Fig 1, P represents 

the input vector with P,, Pz... Pp elements. The vector n is the input vector of each neuron which 

is obtained by multiplication of the input vector elements by the weight of each neuron. The 
k 


weight matrix of the network is w;;, where i represents the number of neurons in each layer, j 


' Tunnel Boring Machine 
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represents the number of input vector elements, and k is the number of layers. Moreover, the 
vector a is the vector of outputs which is calculated by the transfer function f [37]. 


Inputs First Layer Second Layer Third Layer 


al =f1(W'p+b!) wv =f2(Wal+b?) a = £3(Wa2+b!) 
a = £3(W'f2(Wef !(Wip+b!)+b2) +b) 
Fig. 1. Structure of multilayer Artificial Neural Network(ANN) [37]. 


In multilayer feedforward neural networks with backpropagation capabilities, which are 
commonly used in solving rock and soil mechanics problems, neurons are sorted in the form of 
layers including an input layer, an output layer, and one or more middle layer (hidden layers). Fig 
2 shows a schematic view of a neural network. 


4, 

~ ‘ a 
“~ L fl) fl, 

CK hon 


Processing element 


Artificial neural network 
Fig. 2. ANN neurons and calculation process schematic [38]. 


Every processor unit (neuron) in a certain layer is connected to all or some of neurons from other 
layers by weights. The numerical value of the weights controls the intensity of the connection 
between neurons. Zero weight implies no connection among the neurons while negative weight 
implies no connection is allowed among them. Neurons receive weighted inputs from a large 
number of other neurons and send it to the transfer function after adding their bias to it; 
consequently, the output of the transfer function act as the input for the next layer. The bias value 
is applied for correct scale adjustment at neuron input, improving network convergence. 
Considering the artificial neural network in Fig 2, the procedure above could be presented by 
Eqs. 1 and 2 for neuron j [39]: 


Ty = 0; + Vina Wri (1) 


yj = f() (2) 
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in which J; is activation level of neuron j, w;; is connection weighs between j and i neurons, x; is 
inputs from neuron i (i = 0.1.2.....7), 6; is bias of neuron j, y; is output of neuron j, f(.) is 
transform function [37]. 


3. ANN training data generation 


In order to use ANN, firstly it needs to provide a series of data as network input and another 
series of data as network output. In this work regarding of demand data, the common and 
effective methodology was applied in order to create network necessary data. Data generation 
process and steps of data preparation are in followed parts: 


3.1. Discrete fracture network (DFN) generation 


In order to generate the 3D model of the fracture network, the 3D code of 3DFAM was 
developed in MATLAB Programming Language. In 3DFAM Code, the fractures center has been 
defined on the basis of modified Beacher disk and the number of fractures has been generated 
with the homogeneous Poisson process [7]. The uniform distribution has been used in order to 
simulate the fractures location in 3-D [40]. Having generated the 3D block of the discrete 
fracture network, it would be possible to create sections in any direction or even inclined (with 
any dip and dip direction) and any location. The fractures, which are situated outside the study 
limits, are omitted. In addition, the external trace length of the fractures, which some parts of 
them are inside the window, is omitted as well. The main object of truncation of the traces length 
is to leave out the edge effect [3]. 


All the trace lengths existing inside the window under consideration are studied after creating 
sections in the required directions. In order to study the effective connectivity between the 
fractures, it is needed to omit the fractures or cluster of fractures which are hydraulically inactive 
and do not contribute to transferring the flow after creating the section. This process is called 
Fractures Regulation. 


Connectivity is a comprehensive parameter that is a function of orientation, size, spacing, and 
density of the joint sets. Such characteristics are used to define the connectivity portion in the 
rock mass [41]. The flow routes and the number of connected fractures along the flow routes are 
computed by regulating the network. In order to assess the effect of the geometrical 
characteristics of the fractures on the network and the parameters corresponding to the effective 
connectivity as well as the hydraulic conductivity of the fracture network, 3DFAM code can 
compute a number of these parameters for the cases before and after regulating the fractures. The 
percolation frequency (total connectivity) (€), hydraulic conductivity (yn), intersection density 
(P29), intersection intensity (P,,) and termination index (7;) are the most outstanding parameters 
which are going to be calculated after regulating so as to study the status of the fracture network 
effective connectivity. 


The percolation frequency is the very degree of interconnection or total connectivity of a fracture 
network which is computed using Relation 3: 
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= Mnode (3) 


n 


Where Nyoge iS number of intersection nodes and n is total number of fractures in fracture 
network [42]. 


The conductivity parameter, n, represents the connectivity in the fracture network. When a 
uniform distribution characterizes the opening, the conductivity Parameter has a linear 
correlation with the hydraulic conductivity of the fracture networks. This parameter is calculated 
using Relation 4: 


_ nl 1+2nnode/n 
1 = 2 | tnode/m (4) 


Where Nyoge is the number of intersection nodes, n is total number of fractures, A is the model 
area and | is the arithmetic mean of all fracture lengths in fracture network [42]. 


The intersection density (P29) and intersection intensity (P21) represent the connectivity in the 
fracture network. These parameters are equal to the number of intersections and the sum of 
fracture trace length divided by the area of the fracture network. 


3.2. DFN data acquisition 


In order to get hold of the real conditions, four synthetic joint sets are simulated. In 84 
simulations, the input data include Fisher's constant, mean and standard deviation of trace length, 
diameter, aperture and frequency and output parameters include percolation frequency, hydraulic 
conductivity, intersection density (P29), intersection intensity (P2,) and total discharge at outflow 
boundary of fracture network in horizontal (Q,3) and vertical (Q2,) direction flow. 


One of the most important parameters in definition of discontinuities orientation is Fisher's 
constant. K is the coefficient of the Fisher's constant and a positive number which indicates the 
degree of data scattering. The larger amount of this coefficient depicts less scattered data. The 
larger amount of this coefficient depicts less scattered data [43,44]. The Fisher's constant is 
generally in the range of 5 to 440. The range of input data variation are presented in table 1. 


Table 1 
Variation limits of simulations input parameters. 
Fisher Mean Trace Trace Length Standard F Frequenc 
coefficient Length(m) Devistion(m) Dismerer{(ay: Spore pa . 
5 0.77 0.21 0.99 0.0005 0.99 
50 1 0.31 1.25 0.000725 1.25 
100 1.27 0.4 1.6 0.000725 1.59 
250 1.5 0.48 1.86 0.000725 1.86 
250 1.72 0.53 1.93 0.000725 1.93 


In all simulations, the shape of fractures is assumed circular, and DFN is generated in 3-D block 
with dimensions of 5 x 5 x 5 m3. A sample of generated DFN and a 2-D cross-section of the 
model in either non-regulated and regulated form is shown in figure 3. In each block, simulation 
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comprise 3 DFN realizations. Hence in this study that include 84 simulations, 252 realizations of 
DEN is analyzed. 


GER Joint Set: 1 


HE Joint Set: 2 
GE Joint Set: 3 
GGBDBoint set: 4 


on sett ee x (East) 


ed 


if 


ai 
Se 


Fig. 3. ANN Data generation: A. Generated DFN by circular fractures, B. 2-D cross-section of DFN 
(simulated fracture traces) in non-regulated form, C. 2-D cross-section of DFN (simulated fracture traces) 
in a regulated form. 


In order to flow analyzing in DFN, literatures mostly recommended two type of boundary 
condition which is prescribed on the model. In the first type of boundary condition, hydraulic 
gradient is assumed horizontal and flow is calculated based on the outflow of boundary 3. In the 
second type of boundary condition, hydraulic gradient is assumed vertical and flow is calculated 
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based on the outflow of boundary 4. Input and output Pressure head in both condition 
respectively is 5 and | meter. Figure 4 shows above conditions. 
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Fig. 4. Boundary conditions which are applied to DFN floe calculations: A. horizontal hydraulic gradient, 
B. vertical hydraulic gradient [44]. 


4. Optimum ANN design 


In optimal neural network designing, the correct selection of input variables has a large impact 
on the efficiency of the network [36]. In this study, six parameters (Fisher coefficient, the length 
of the fractures, the standard deviation of the length of the fractures, the diameter of the fracture, 
the number of discontinuities in the unit length of the scan line and the opening of the fractures) 
were assumed as the inputs of the neural model. Given the number of input parameters, the input 
layer has six neurons. On the other hand, the outputs of the neural model have six parameters, 
namely: density of fracture intersections (P29), the intensity of fracture intersection (P21), 
frequency of permeation (C), conduction parameter (CP), the average exiting flow rate on the 
boundaries normal to the gradient direction and parallel to the horizontal flow (Q,3) and parallel 
to the vertical flow (Q2,4), that are hydraulic parameters of the discrete fracture network. Given 
the number of network outputs, the output has 6 neurons. The neural network that was employed 
in this study is a multilayer feedforward network with backpropagation. Using the neural 
network toolbox of MATLAB, 50 various models of neural networks with different structures 
were investigated. Ultimately, the neural network model in Fig. 5 yielded the best responses with 
the least error and the highest correlation between model outputs and target values. Parameters 
and characteristics of the optimal neural network is presented in Table 2. 


Hidden Layer Output Layer 


10 6 
Fig. 5. Optimal designed ANN blocks. 
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Gradient Moments decreasing can train any network as long as its weight, net input, and transfer 
functions have derivative functions. Backpropagation is used to calculate derivatives of 
performance perf concerning the weight and bias variables X. Each variable is adjusted 
according to gradient descent with momentum: 


dX = mc x dXprev + Ir *(1— mc) * dperf /dxX (5) 


where dXprev is the previous change to the weight or bias, parameter Ir indicates the learning 
rate, derivatives of performance dperf with respect to the weight and bias variables X and 
parameter mc is the momentum constant that defines the amount of momentum. mc is set 
between 0 (no momentum) and values close to 1 (lots of momentum). A momentum constant of 1 
results in a network that is completely insensitive to the local gradient and, therefore, does not 
learn properly. 


Training stops when any of these conditions occurs: 
*The maximum number of epochs (repetitions) is reached. 


*The maximum amount of time is exceeded. 


*Performance is minimized to the goal. 


Table 2 
Characteristics of optimal designed ANN. 


Levenberg — Marquardt Learning Function 


Gradient Moments decreasing Weighs and bias training Algorithm 


78 Total data number 

2 Network layer Number 

6 Input Layer Neurons Number 
10 Hidden Layer Neurons Number 
6 Output Layer Neurons Number 


5. Results and discussion 


Various methods can be used to evaluate the performance of the neural network. Mean squared 
error diagrams, error histograms, and also correlation coefficient between the outputs of the 
neural network and actual values were used in this study for evaluation of the designed neural 
network performance. The neural network performance was evaluated through investigation of 
mean squared errors by common techniques. Fig 6 shows the mean squared errors calculated for 
every frequency of the alternative neural network computation along with its descending 
development for the optimal designed neural network. 
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2 Best Validation Performance is 2.287 at epoch 7 


Mean Squared Error (mse) 


13 Epochs 


Fig. 6. Optimal designed ANN calculation epochs Mean Squared Error(MSE). 


Experimental and evaluation data curves showed similar behavior proving the stability of the 
designed network and its reliability. According to Fig 6, the best (the least) mean squared error of 
evaluation data, being 2.278, was obtained at the 7th stage of the learning process. If the test data 
curve showed a significant increase ahead of the evaluation data, it would have implied a fault in 
network performance which is not seen in Fig 6. 


Fig 7, shows the difference histogram for actual and computed values. Blue bars in the figure 
show learning data error, green bars show evaluation data error and the red bars show the 
network test data error. The error histogram shows the majority of error values to be between -3 
and 2. Fig 7 shows that more than 90% of the data have very low errors in the range of -1.16 to 
0.85. Only a few data are associated with a larger error, showing the accuracy of performance. 


Error Histogram with 20 Bins 


EEE Training 
420 HE validation 
EE Test 


Zero Error 


Instances 


-19.19 
17.18 
-15.18 
13.18 
14.17 
9.169 
-7.166 
5.163 
3.159 
1.156 
10.86 
12.87 
14.87 
16.87 
18.88 


0.8474 
2.851 
4.854 
6.858 
8.861 


Errors = Targets - Outputs 


Fig. 7. Optimal designed ANN Error Histogram. 
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The most important metric for evaluating the performance of the neural network is the 
correlation between actual target data and data computed by the neural network. Fig 8 shows the 
correlation between actual target data and computed data for all three groups of learning, 
evaluation and test data. 
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Fig. 8. Optimal designed ANN Error Histogram. Correlation Coefficient between predicted and actual 
data. 


Diagrams in Fig 8 show an excellent correlation coefficient between actual target data and 
computed data. This shows the desirable performance of the designed neural network. 
Correlation coefficient values are very close to 1 and are presented in Table 3. 
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Tabl 
ie ee between actual and predicted value using optimal designed ANN. 
The relation between predicted and actual data Correlation Coefficient Data type 
Pridicted value = 0.99 * Target + 0.6 0.99 Training data set 
Pridicted value = 0.99 « Target + 0.46 0.99 Validation data set 
Pridicted value = 0.97 * Target + 0.48 0.97 Test data set 
Pridicted value = 0.99 * Target + 0.62 0.99 All data 


6. Conclusion 


In this paper, Given the importance of hydraulic parameters of the discrete fracture network, and 
their time-consuming computation for dense real models, a predicting model was designed and 
developed drawing on the capabilities of the artificial neural network. The model predicts six 
hydraulic parameters upon receiving six input parameters, namely: The Fisher coefficient, length 
of the fractures, the standard deviation of the length of the fractures, the diameter of the fracture, 
the number of discontinuities in a unit length of the scan line and opening of the fractures. The 
hydraulic parameters are: the density of fracture intersection (P20), the intensity of fracture 
intersection (P21), frequency of permeation (¢), conduction parameter (CP), average exiting flow 
rate on boundaries normal to fluid gradient direction (Q13), and the average existing flow rate on 
boundaries parallel to the fluid gradient (Q24). The discrete fracture network model was created 
using the developed 3DFAM software with the help of which neural network output and input 
parameters were computed. The high accuracy of the designed neural network provides a more 
accurate analysis of the hydraulic behavior of the discrete fracture network by predicting the 
aforementioned hydraulic parameters. Robust fitting tools that are capable of receiving multi- 
dimensional inputs and generating multi-dimensional outputs is a characteristic of this network. 
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